function DJ = s_grad(model,fid,kbasis,cgp)

% S_GRAD Compute gradient of error function
%
% MODEL is the parameter array
% FID is the FID array
% KBASIS is a basis array defined by make_kbasis().
% CGP is the reconstruction parameter structure.

if isfield(fid,'N') 
    NF = fid.N;
else
    NF = length(fid.sig);
end
if isfield(model,'N') 
    NM = model.N; 
else
    NM = length(model.pvec)/2; 
end
sigdiff = fid.sig(1:NF) - fid.est(1:NF);
%figure(16)
%plot(real(sigdiff))
%drawnow
wfr = model.wfr;
wfi = model.wfi;
ix = model.ix;
iy = model.iy;
alphaM = cgp.alphaM;
alphaf = cgp.alphaf;
if (isreal(alphaf))
    alphaf = (1+1i)*alphaf;
end

p_image = zeros(model.reso);
p_grad = zeros(model.reso);
p_image(ix+model.reso*(iy-1)) = model.pvec(1:NM);
M_grad = zeros(NM,1);
f_grad = zeros(NM,1);
CC_conj = conj(-sigdiff);
L = [ 0 1 0; 1 -4 1; 0 1 0];
M = transpose(model.pvec(1:NM));

wf = weight_cexp(model.pvec(NM+1:2*NM).',wfr+1i*wfi);
wf_n = ones(NM,1);
if isstruct(kbasis)
    % use separable calculation
    M_grad = (kbasis.ekxx(1,ix).*kbasis.ekyy(1,iy)).'*CC_conj(1);
else
    M_grad = kbasis(:,1)*CC_conj(1);
end

if (alphaM ~= 0)
    p_grad = conv2(conv2(conj(p_image),L,'valid'),L);
    M_grad = M_grad + alphaM*p_grad(ix+model.reso*(iy-1));
    if (1==0)
        figure(16)
        subplot(1,2,1)
        imagesc(real(p_image))
        subplot(1,2,2)
        imagesc(real(p_grad))
        drawnow
        pause;
    end
end

for ii = 2:NF
    wf_n = wf_n.*wf;
    if isstruct(kbasis)
        % use separable calculation
        Cekxx = kbasis.ekxx(ii,:)*CC_conj(ii);
        temp = wf_n.*(Cekxx(ix).*kbasis.ekyy(ii,iy)).';
    else
        temp = wf_n.*kbasis(:,ii)*CC_conj(ii);
    end
    M_grad = M_grad + temp;
    f_grad = f_grad +(ii-1).*temp;
end

f_grad = M.*f_grad;
if (alphaf ~= 0)
    p_image(ix+model.reso*(iy-1)) = model.pvec(NM+1:2*NM);
    p_grad = conv2(conv2(conj(p_image),L,'valid'),L);
    p_grad = wfr*real(p_grad) + wfi*imag(p_grad);
    pg = p_grad(ix+model.reso*(iy-1));
    f_grad = f_grad + wfr*real(alphaf)*real(pg)+1i*wfi*imag(alphaf)*imag(pg);
end

DJ = [2*conj(M_grad); 2*wfr*real(f_grad)-1i*2*wfi*imag(f_grad)].';
